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ABSTRACT 

We propose a closure model for the transport of entropy and momentum in astro- 
physical turbulence, intended for application to rotating stellar convective regions. 
Our closure model is first presented in the Boussinesq formalism, and compared with 
laboratory and numerical experimental results on Rayleigh-Benard convection and 
Homogeneous Rayleigh-Benard convection. The predicted angular momentum trans- 
port properties of the turbulence in the slowly rotating case recover the well-known 
A— effect, with an amplitude uniquely related to the convective heat flux. The model 
is then extended to the anelastic case as well as the fully compressible case. In the 
special case of spherical symmetry, the predicted radial heat flux is equivalent to that 
of mixing-length theory. For rotating stars, our model describes the coupled transport 
of heat and angular momentum, and provides a unified formalism in which to study 
both differential rotation and thermal inhomogeneities in stellar convection zones. 

Key words: convection; hydrodynamics; turbulence; stars:rotation 



1 INTRODUCTION 

Turbulent convection occurs frequently in stellar interiors 
and other astrophysical fluid flows. While convective mo- 
tion naturally transports heat and chemical elements, the 
transport of angular momentum by convection in rotating 
bodies is a more subtle issue. It is of particular interest in 
the case of the Sun, where the internal pattern of rotation 
has been measured but remains incompletely understood. It 
may also play an significant role in accretion flows. 

Numerical simulations of astrophysical convection are 
becoming increasingly powerful and capable of resolving a 
widening range of length and time scales. Nevertheless, a 
simpler, statistical description of turbulent transport is de- 
sirable in order to treat the effects of convection on the struc- 
ture and evolution of stars. It almost goes without saying 
that such a description cannot be derived strictly from the 
equations of fluid dynamics but must involve some modelling 
or parametrization. 

The mixing-length theory of turbulent transport was 
developed by Prandtl (1925) and applied to stellar convec- 
tion by Biermann (1932). It is still the basic model used in 
most calculations of stellar structure and evolution, usually 
in the form devised by Bohm-Vitense (1958). The main pur- 
pose of mixing-length theory is to relate the convective heat 



flux to the superadiabatic gradient; in this context it does 
not usually deal with the transport of (angular) momentum 
that arises in the presence of shear or rotation. 

A standard theoretical approach to convection in dif- 
ferentially rotating stars is set out in the monograph by 
Riidiger (1989). Angular momentum transport is described 
by a Reynolds stress tensor whose components can be re- 
lated to the large-scale mean flows and thermodynamical 
gradients. A first contribution to the Reynolds stress is typ- 
ically proportional to the angular velocity gradient through 
a turbulent viscosity coefficient. An important additional 
contribution comes from the A-effect (named after Lebedin- 
sky), whereby even uniformly rotating convection transports 
angular momentum by virtue of its anisotropy. Attempts to 
constrain or parameterize these quantities have been made 
through local numerical simulations (e.g. Kapyla, Korpi & 
Tuominen 2004) or theoretical models (e.g. Kitchatinov & 
Riidiger 1993). Mean-field models of stellar rotation (e.g. 
Kitchatinov & Riidiger 1999, Rempel 2005) have been de- 
veloped which use such parameterized expressions for the 
Reynolds stress and heat flux. 

Reynolds-stress models of turbulent flows have been de- 
veloped in the engineering community over several decades 
(e.g. Pope 2000). The exact equation governing the Reynolds 
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stress in a turbulent fluid cannot be solved because of the 
well known closure problem whereby an infinite hierarchy 
of correlations is involved. Nevertheless, by parametrizing 
the difficult terms in this equation, models can be con- 
structed that bear some fidelity to the turbulent dynam- 
ics. Prom a more physical point of view, what is obtained 
is a time-dependent constitutive equation for the turbulent 
fluid, which relates the turbulent stress to the local history of 
deformation. There is a close similarity with models of non- 
Newtonian fluids (Ogilvie & Proctor 2003). The advection 
and deformation of the turbulent stress are accurately repre- 
sented since they derive from linear terms in the Reynolds- 
stress equation, while the nonlinear 'relaxation' effects are 
only modelled (as is also true for non-Newtonian fluids). 

A similar approach can be applied to turbulent convec- 
tion in which buoyancy forces play an essential role. The ad- 
ditional correlations that must be considered are the flux and 
the variance of entropy (or temperature, in the Boussinesq 
approximation). This approach offers some benefits over the 
conventional description in terms of a turbulent viscosity 
and a A-effect. It can be formulated in a covariant manner 
and is not tied to the spherical geometry of a slowly rotat- 
ing star. It starts from a more fundamental description and 
allows phenomena such as the A— effect to emerge in a nat- 
ural way from more elementary considerations. It may also 
allow a more unified approach to be taken towards problems 
involving astrophysical turbulence. 

In this paper we explore some of the consequences 
of a simple dynamical model of astrophysical convection 
of this type. The model derives from one originally con- 
ceived for magnetohydrodynamic turbulence in accretion 
discs (Ogilvie 2003) and later applied to rotating shear fiows 
without magnetic fields (Garaud & Ogilvie 2005, GO05 her- 
after). Our motivation is to develop and test a model that 
can be applied to the convcctive zone of the Sun, to other 
stars or to accretion discs. We emphasize, however, that our 
model is chosen to be as simple as possible for the purposes 
of this investigation. In contrast with some of the engineer- 
ing literature, we restrict the algebraic complexity in order 
to retain a physical understanding of the terms in the equa- 
tions. Purther refinements are likely to be required in order 
to provide an accurate match to a wide range of data. 

In comparing a closure model of astrophysical convec- 
tion with experimental and numerical results, we face certain 
difficulties. Astrophysical convection usually takes place at 
very high Rayleigh number, in a highly turbulent regime. 
Experiments have been conducted at very high Rayleigh 
number but mainly for the Raylcigh-Bcnard problem in 
which the fiow is dominated by boundary layers, which may 
not be relevant in an astrophysical context, or by mean fiows 
not represented in the closure model. An alternative system 
is provided by the homogeneous Rayleigh-Benard problem, 
which has periodic boundary conditions in all directions. 
This model, however, has certain peculiarities of its own. 
These issues will be addressed in the sections that follow. 

In the remainder of the paper, we develop the closure 
model first in the Boussinesq approximation (Section 2) 
and apply it to the standard Raylcigh-Bonard problem 
(Section 3). We then consider the homogeneous Rayleigh- 
Benard system with triply periodic boundary conditions 
(Section 4); in this section we also introduce rotation and 
discuss the A— effect. We then adapt the model to the anelas- 



tic approximation for use in stars and other astrophysical 
fiows (Section 5) and finally draw conclusions (Section 6). A 
number of technical details are covered in the appendices. 



2 CLOSURE MODEL IN THE BOUSSINESQ 
SYSTEM 

2.1 Basic equations 

In the Boussinesq approximation (e.g. Chandrasekhar 1961) 
the equations governing the motion of the fluid are 



diUi = 0, (1) 

po{dt + Ujdj)ui = pQi - dip + povdjjUi, (2) 

p = po[l-a(T-To)], (3) 

{dt + Uidi)T = KdiiT, (4) 



where we have adopted a Cartesian tensor notation. The 
dynamical variables are the velocity u, the density p, the 
pressure p and the temperature T. Quantities regarded as 
constant in the Boussinesq approximation are the reference 
density po, the reference temperature To, the coefficient of 
expansion q, the gravitational acceleration g, the kinematic 
viscosity v, and the thermal diffusivity k. 

A simple, static basic state is possible when the temper- 
ature is uniform and the pressure gradient balances gravity, 
i.e. 

T = To, (5) 

p = Po + PogiXi, (6) 

where po is a reference pressure. To examine departures from 
this state we define 

e = T - To, (7) 
po 

obtaining the governing equations 

diUi = 0, (9) 

{dt + Ujdj)ui = -a&gi - diip + vdjjUi, (10) 

{dt + Uidi)Q = ndii®. (11) 

2.2 Fluctuations 

Wc now adopt a standard procedure and separate the dy- 
namical variables into mean and fiuctuating parts, e.g. 

Ui = Ui +u'i, {u'i) = 0, (12) 

where the angle brackets or the overbar are interchangeably 
used to denote a suitable averaging operation such as a tem- 
poral, spatial or ensemble average. The mean parts of the 
governing equations are 

diUi = 0, (13) 
{dt + Ujdj)ui = —aQgi — di%l) -\- vdjjUi — djRij, (14) 
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{dt + Uidi)e = KdiiQ - diFi, (15) 

where 

Rij = u'iu'j (16) 

is the Reynolds tensor, representing (minus) the turbulent 

stress, and 

F^ = e'u'i (17) 

represents the turbulent heat flux density. The problem at 
hand is to determine Rij and Fj and thereby close the system 
of mean equations. We also introduce the quantity 

Q = e'^ (18) 

representing the temperature variance. It should be noted 
that all three quadratic correlations Rij, Fi and Q will be 
redefined when we move on to the (more relevant) anelastic 
system in which the reference density is non-uniform, but 
these definitions are convenient for the Boussinesq system. 
The fluctuating parts of the governing equations are 

diu'i = 0, (19) 

(dt + Ujdj)u'i + u'jdjUi = —aQ'gi — di-ip' + vdjju'i 

-dj{Rij - Rij), (20) 

{dt + Uidi)e' + u'idiQ = KdiiQ' - di{F, - Fi). (21) 

Prom these we can obtain exact equations for Rij, Fi and Q 
in the form 

{dt + Ukdk)Rij + RikdkUj + RjkdkUi 

+a{Figj + FjQi) - vdkkRi] = -{u'idjtp' + UjO^-ij/} 
-{u'idkRjk + u'jdkRik) - 2i^{dku[dku'j), (22) 

{dt + Ujdj)Fi + RijdjO + FjdjUi + aQgi - ^{1/ + K)djjFi 
= ~{e'di^p') - {Q'djRij + u'idjFi) 
+\{u-K){dj{Q'dju'i-u'idjQ')) 
-{v + K){dju'idje'), (23) 

{dt + Uidi)Q + 2FidiQ - KdiiQ 

= -2{e'diFi)-2K{{diQ'f). (24) 

The left-hand sides of these equations represent the linear 
interaction of Rij, Fi and Q with the mean velocity gradi- 
ent, the mean tcmijerature gradient and the gravitational 
field, as well as their diffusion by the microscopic transport 
coefficients. There is no difficulty in treating such terms ex- 
actly as they appear. The right-hand sides of these equations 
contain difficult terms of three sorts: those involving correla- 
tions with the pressure fiuctuation ip' , those involving triple 
correlations of fiuctuating quantities, and dissipative terms 
involving the microscopic diffusivities v and k. These effects 
can all be regarded as 'non-linear'; although viscous diffu- 
sion, for example, is a linear process, when the Reynolds 
number is large the viscous terms can be significant only 
when a turbulent cascade has forced structure to appear on 
the dissipative scales. None of the terms on the right-hand 
sides of these equations can be written in terms of Rij, Fi 
and Q without further knowledge of the statistical proper- 
ties of the fluctuating quantities, such as the spectrum of the 
turbulence, which are determined by the non-linear physics 
of the turbulent cascade. 



2.3 Proposed closure model 

We therefore attempt to model the system by retaining the 
exact forms of the left-hand sides and proposing simple clo- 
sures for the right-hand sides, i.e. 

{dt + Ukdk)Rij + RikdkUj + RjkdkUi 
+a{Figj + FjQi) - udkkRij 

= J-ij {Rij ,Fi,Q, . . .), (25) 

{dt + Ujdj)Fi + RijdjQ + FjdjUi + aQgi - 5(1^ + K)djjFi 
= J^i{Rij,Fi,Q,...), (26) 

{dt + Uidi)Q + 2FidiQ - Kdi.Q = T{R,j,Fi, Q,...), (27) 

where the quantities T are non-linear tensorial functions of 
their arguments. The dots represent the parameters of the 
problem, on which the functions T may depend. 
A simple example of such a model is 

{dt + Ukdk)Rij + RikdkUj + RjkdkUi 
+a{Figj + Fjgi) - vdkkRij 

= — —R}^^Rij — —R}^^{Rij — ^R5ij) — v^^Rij, 

(28) 

{dt -\- Ujdj)Fi + Rijdj® + FjdjUi + aQgi - ^{1/ + K)djjFi 
= _^Ry^Fi-\{v + n)^F, (29) 

{dt+Uidi)Q + 2Fidi@- KdiiQ = ~R^''^Q-K^Q, (30) 

where R = Ra is the trace of the Reynolds tensor, which 
is twice the turbulent kinetic energy per unit mass, and Ci, 
C2, Cq and C7 are positive dimensionless coefficients of order 
unity, of a universal nature. (Coefficients C3, C4 and C5 
are reserved for a magnetohydrodynamic extension of the 
model, see Ogilvie 2003) 

The justification for introducing non-linear terms of the 
above form is similar to that used in the model of mag- 
netorotational turbulent stresses originally introduced by 
Ogilvie (2003). The term involving Ci causes a dissipation 
of turbulent kinetic energy, and allows for the free decay 
of hydrodynamic turbulence. The term involving C2 redis- 
tributes energy among the components of Rij, and corre- 
sponds to the tendency of hydrodynamic turbulence to re- 
turn to isotropy through the effect of the pressure-strain 
correlation. Both arc constructed assuming that these ef- 
fects occur on a tiinescalc related to the eddy turnover time, 
L/R}^^ , where L is defined as the typical scale of the largest 
turbulent eddies. Terms Ce and C7, related to the transport 
of heat, are advanced by simple analogy. The coefficients 
must satisfy certain conditions to ensure the rcalizability of 
the model, as discussed in Appendix A. 

The terms proportional to the microscopic diffusion co- 
efficients are introduced to allow a modelling of the correla- 
tion terms 2v {dku'idku'j) , {v + K){dju'idjQ') and 2K{{diQ')'^) 
at moderate Reynolds number, i.e. close to the onset of con- 
vection. In such a situation a turbulent cascade does not 
form and the dissipative terms are proportional to, rather 
than independent of, the diffusion coefficients. In a similar 



© 0000 RAS, MNRAS 000, 000-000 



4 P. Garaud, G. I. Ogilvie, N. Miller, S. Stellmach 



way, for turbulent shear flows, GO05 proposed to model the 
momentum diffusion term as 



(31) 



on dimensional grounds. Indeed, it is expected that near the 
onset of convection, most fluid motions will be on the largest 
scales of the system (L) . By analogy, we model the other two 
terms here as 



(32) 
(33) 



{ly + K){dju'idjQ') - 

Therefore the dissipative term in each of equations (I25|l - (|27|) 
is modelled by a sum of two terms, one that is independent 
of the diffusivity and dominates at high Reynolds numbers, 
and another that is proportional to the diffusivity and dom- 
inates at moderate Reynolds numbers. This completes the 
justification for the form of the closure model proposed in 
equations ([28]), (l29j and (|30)) . 



where g = —g^. In the case of no-slip boundaries with fixed 
temperature on each plate as listed above, R, Rzz, and 
Q are zero on both boundaries. 

This system of ODEs with associated boundary condi- 
tions can be solved with a two-point boundary-value solver. 
Typical solutions are shown in Fig. [1] for various Rayleigh 
numbers, defined here as 



Ra : 



agh^AT 



We set the Prandtl number 



Pr : 



(35) 



(36) 



to 1 for the purposes of illustration. Note the appearance of 
the characteristically fiat temperature profile between the 
two plates as Ra — >■ oo and of the thin thermal boundary 
layers. We now study in more detail the structure of the 
solution. 



3 RAYLEIGH-BENARD CONVECTION 
3.1 Model setup 

We now apply the closure model to the problem of Rayleigh- 
Benard convection. We consider a horizontally infinite, 
plane-parallel system, where the bottom plate is located at 
height z — and the top plate at height z — h. The relative 
temperature of the bottom plate is O = AT while that of 
the top plate is O = 0. 

In this setup, we look for statistically steady and hori- 
zontally homogeneous solutions assuming that mean quan- 
tities and correlations between fluctuating quantities vary 
only with z. We also assume that there are no mean flows 
in the system. Equations (IT51l - (fTS)l and pSjl - PUIl reduce to 
a set of ordinary differential equations (ODEs) which can 
be solved to obtain the temperature profile 0{z) between 
the two plates, the profiles of the turbulent kinetic energy, 
R{z)/2, and the temperature variance, Q{z) (for example). 

By analogy with Prandtl's mixing-length formulation 
(Prandtl, 1932) we set L, the size of the largest eddies, to 
be equal to the distance to the nearest wall, i.e. L{z) = 
min(2:, h—z) (see GO05 for applications of the same principle 
to pipe fiows and to Couette-Taylor fiows). 

It can be shown with little effort that R^y = R^z = 



Rv 



0, as well as 



0. The remaining set of five 



second-order ODEs fully characterizes the system: 



,'13. = ^^r+^r:^/^ 



dz^ 
A^R.z 



L2 

-2aFzg, 



L 



2aFzg, 



dz^ 



d^Q 

^dz^ 

d^e 
'di^ 



l(^ + '^) ^2 

-aQg + Rz2 



de 

dz' 



dFz 
dz ' 



+ 2F, 



de 

dz' 



(34) 



3.2 Universal profile of convection from a wall 

As in the case of shear flows past a wall (see GO05), we can 
derive a universal proflle for convection away from a wall. 
Let us consider a semi-inflnite domain z > 0, in which case 
L = z, and let Fq be the convective heat flux through the 
system. We define dimensionless variables via 

,2., 1 1/4 



_ agFo 
F. = Fofiv), 



Q 



6- AT: 



F^v 



-I 1/4 



agK,^ 

iv), 



1/2 



liv), 



R^ 



agFoK 



21 1/2 



(37) 



so that the system of equations (|34p becomes 

// 1 Cl 3/2 „/• 

r = —r + r ' - 2f, 



Pr + 1 



f" = 



n 
a 
v' 

Pr + 1 C, 



,2 + Pj. 



Pr 1] 

1 Cl + C2 1/2 

r ' r 



1 C: 



Pr 3v 



2 3/2 
-r ' 



2/, 



Ck I 1/2 I r,rn' 

— gH r q + 2fe , 

^2 



(38) 



The boundary conditions at ?7 = are r — Tzz ~ f = q = 
9 = 0. 

Solutions very close to the wall {rj <^ 1) satisfy: 



r and rzz oc rj°'" with a,^{av — 1) = C'v, 
f oc rj°""^ with a^f^ia^K ~ 1) = C^«;, 
q oc r;"" with ««(«« — 1) = Ck. 



(39) 



These simple relationships provide an ideal way of calibrat- 
ing each of the three constants C,^, C^^ and Ck individually 
(see Section [3.4|) . by analysing the power-law behaviour of 
the near-wall profiles of experimental or numerical data. 
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Figure 1. Vertical profiles @{z) (in units of AT), R{z) (in units 
oi K^/h^), F:,(z) (in units of kAT/Zi) and Q(^) (in units of (AT)^) 
for Ra = 10*^ (dotted line), Ra = 10* (dashed line) and Ra = 10^" 
(solid line). In all cases, Pr = 1. 



Solutions far away from the boundary layer can be ex- 
panded as 



e = eo + 3/ir,-^/=' + 0(r?-^/^), 



where 
ro = 
/i = 



(f) 



2/3 



-1/2 



3Ci + C2 



3(Ci 



2/i 



'"0, 



Ci 



1/2 ' 



(41) 



C7 ^ 3(crTcJT 

However, unlike ro, /i and go the constant 6*0 cannot be de- 
termined without a numerical calculation of the boundary- 
layer solution for 77 = 0(1). 

The scaling laws obtained for r, f, 6 and q far from the 
wall are expected on dimensional grounds, and recover the 
well-known solution of Priestley (1954). They are analogous 
to the universal "log-law" solutions for turbulent shear flows 
past a wall (e.g. Schlichting, 1979). By comparing profiles of 
r, / and q with laboratory or numerical experiments, one 
can constrain some of the unknown coefficients {Ci} (see 
Section Ea. 



3.3 Nusselt— Rayleigh number relationship 

The heat flux through the system in Rayleigh-Benard con- 
vection is commonly measured by the dimensionless Nusselt 
number 



Nu = l-f 



hFo 
kAT' 



(42) 



which compares the total heat flux with the conductive one 
in the absence of convection. The universal convection-from- 
a-wall solution calculated in the previous section can be used 
to derive the relationship between the Nusselt number and 
the Rayleigh number. 

Indeed, by selecting a Rayleigh number we set the rela- 
tive temperature at the midpoint z — h/2 to be 6 = AT/2 
which implies through (|37|) that 



(43) 



yielding an equation for the (unknown) constant heat flux 
Fq. In dimensionless terms, we have the implicit equation 
for Nu: 



AT 


- AT = 


r 1 


1/4 / 

• 


agFo 


1/4 















-[Ra(Nu- 1 



,-3il/4 



(^|[Ra(Nu- 



11/4 



(44) 



which can be solved to find Nu(Ra). In the limit of very 
large Rayleigh number the mid-point of the system is very 
far from the boundary layer, so ~ 60 which then recovers 
the standard scaling law (Malkus 1954) 



Nu = 1 -f 



/ Ra 



1/3 



(45) 



(40) 



The constant 60 depends only on Pr and on the closure pa- 
rameters {Ci}, but cannot easily be expressed analytically 
in terms of these parameters. 



3.4 Comparison with data and estimation of the 
model parameters 

The aim of this section is to estimate, in a rough sense, the 
parameters {Ci} by comparing the model predictions with 
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numerical simulations and laboratory experiments. This ap- 
proach was successfully used in GO05 on pipe flow data and 
Couette-Taylor data, yielding: 

Ci ~ 0.4 , C2 ~ 0.6 , ~ 12. (46) 

Under the assumption that the closure parameters are uni- 
versal properties of the turbulent cascade, these estimated 
values should also apply to the case of turbulent convection 
without need for re-calibration. The remaining parameters 
Ce, C7, C'vK and Ck may then be independently estimated. 
In the following sections, we first discuss this assumption in 
the light of known model limitations. We then select appro- 
priate experimental datasets and use them to constrain the 
remaining parameters. 

3.4-1 Discussion of the model limitations 

As discussed by Ogilvie (2003) and GO05 the closure model 
proposed has two intrinsic limitations: it ignores some 
(but not all) of the effects of pressure-strain correlations 
< u'idjijj' >, and assumes that the effect of all modelled 
terms (such as the triple-correlations in (|28|) - (|30| )) is local 
both in time and space. As a result, it may poorly represent 
strongly sheared systems or systems where the turbulent ed- 
dies exhibit a strong degree of spatial or temporal coherence. 

The neglected effects of the pressure-strain correlations 
are not thought to be important in turbulent convection, 
except in the presence of strong rotation or of an externally 
driven strong mean shear (where the timescale of rotation 
and shear is comparable to that of the convection). The 
closure should be well-suited to model convection in stellar 
interiors, but maybe less so for convectively unstable accre- 
tion discs. We defer this particular case to subsequent work. 
However, for similar reasons these effects are also likely to 
be important in pipe flows or Couette-Taylor flow, which 
were used as a basis for calibrating the constants Ci and 
C2 (see GO05). Consequently, the estimates given in ((46)) 
could be somewhat biased, in particular C2 which contains 
information on the rate of return to isotropy. Comparing 
the model with turbulent convection experiments (see be- 
low) can therefore help refine the estimates for Ci and C2 
using more appropriate data. 

As mentioned above, the closure is also less reliable 
when applied to systems where the turbulence exhibits co- 
herence over large scales or long timescales. This might 
pose some problems when applied to convection in a finite 
domain, since large-scale coherent plumes which span the 
whole system are commonly observed in most cases ranging 
from Boussinesq to fully compressible systems. Comparisons 
with experiments can help reveal which aspects of convec- 
tive transport are adequately described by the model, and 
which are not. 

3.4-2 Available experimental data 

Our application of the closure model to Rayleigh-Benard 
convection in Sections 13.21 and 13.31 assumes for simplicity 
that the system is horizontally invariant, while all labora- 
tory and numerical experiments have a limited horizontal 
extent. The presence, nature and geometry of the side-walls 
are known to affect various properties of the turbulent con- 
vection, in particular through the generation of large-scale 



circulations (often called "wind"). This wind influences the 
overall heat transport properties by changing the nature of 
the boundary layers (Castaing et al. 1989; Cioni et al. 1997; 
Grossmann & Lohse 2000, 2001, 2002, 2004). It also induces 
large-scale horizontal inhomogeneities, so that the measured 
vertical profiles of mean quantities and higher-order mo- 
ments may vary with position (Maystrenko, Resagk & Thess 
2007). While our formalism can in principle be applied to 
finite geometries and self-consistently model the effect of 
large-scale fiows, such an extension is beyond the scope of 
the present paper. 

In order to minimize the effect of side-walls we restrict 
the model comparison to experimental setups with very large 
aspect ratios (defined as the ratio of the horizontal to verti- 
cal extent of the domain, and denoted as F). There are a few 
large aspect ratio, high Rayleigh number experimental stud- 
ies which provide measurements of the Nusselt number. Of 
particular interest are results of Fiinfschilling et al. (2005) 
for convection in water (Pr = 4.38) in a cylindrical enclosure 
of aspect ratio up to F = 6, for Ra up to a few times 10^". 
Niemela & Sreenivasan (2006) provide similar information 
for convection in Helium (0.7 < Pr < 8) in a cylindrical con- 
tainer with F = 4, for Rayleigh numbers between 10* and 
10^^. Finally, the Ilmenau barrel experiments of DuPuits, 
Resagk & Thess (2007) provide Nu(Ra) for convection in 
air (Pr = 0.7) in a cylindrical enclosure with variable aspect 
ratio up to 11.3, for Rayleigh numbers up to a few times 10* 
(in the case of the largest aspect ratio). 

By contrast, only very few large aspect ratio experimen- 
tal measurements of the boundary-layer profiles of velocity 
and temperature correlations (such as Rij, Fi or Q) have 
been reported. The largest aspect ratio experiments avail- 
able (F = 11.3) with fully resolved boundary layer profiles 
are presented by DuPuits, Resagk & Thess (2007) although 
the data provided is limited to the mean and rms tempera- 
ture profiles. 

Taking a different approach, direct numerical experi- 
ments are a powerful tool for "idealized" experiments. Hor- 
izontally periodic simulations minimize the effect of side- 
walls (although retain a finite aspect ratio) and permit re- 
solved and precise measurements of all desired mean and 
fiuctuating quantities within the fiow. The main drawback 
is the limited range of parameter space for which resolved 
simulations can be run (typically, Ra < 10* for large aspect 
ratio simulations at Pr = 0(1)). 

For these reasons, we use a combination of experimental 
data (DuPuits, Resagk & Thess 2007) and numerical sim- 
ulations to calibrate the remaining model parameters. Our 
numerical simulations are all run for Pr = 1, in a horizon- 
tally periodic domain with aspect ratio Lx/L^ = LyjL^ = 4, 
using a spectral method briefly described in Appendix B. 
The largest Rayleigh number achieved in this case is Ra 
= 2.1 X 10^ Fi gure [5] shows a typical snapshot of the re- 
sults, in this parameter regime, for the temperature fleld for 
example. The results of the simulations are globally consis- 
tent with those of Hartlep (PhD thesis, 2005, Gottingen). 

3.4.3 Near-wall profiles and estimation ofCu, Ck, and 

a. 

Very close to the wall (?? ^ 1), the closure model solutions 
for the normalized correlations r, r^^, f, q and 6 are well 



© 0000 RAS, MNRAS 000, 000-000 



Turbulent convection 7 




Figure 2. Volume-rendered visualization of the temperature field 
in our numerical simulation of Rayleigh-Benard convection for 
Ra = 2.1 X 10^, and Pr = 1. The system is doubly-periodic in the 
horizontal direction, with aspect ratio 4, and has no-slip bound- 
ary conditions at the top and bottom boundary. The colour and 
opacity scheme has been selected to emphasize structures near 
the lower boundary layer. 

approximated by power laws, as described in equation (|39p . 
These relationships can be compared with data and provide 
a simple way of individually estimating each of the model 
constants C^, Ck and Cvk from laboratory or numerical ex- 
periments. 

Comparisons of 1)39^ with the experimental near-wall 
profile for rzz{ri), f{ri) and q{ri) yield slopes close to 4 
(see Fig. [3|, close to 3 (see Fig. |4|), and Qk close to 2 (see 
Fig. E]). Note that while the amplitude of the power-law ob- 
served in the near-wall profile for q{ri) is seen to depend on 
the experiment considered, the slope Ok appears to be uni- 
versal. We then adopt the following values for the constants 
Cf, and C^'- 

a = 12 ± 1 , 
= 6 ± 0.5 , 

C« = 2 ± 0.2 . (47) 

Given the experimental and model uncertainties, these val- 
ues and their errorbars should be thought of as rough esti- 
mates rather than precise calibrations. 

It is comforting to note that this independent compar- 
ison recovers the value of found by GO05. Moreover, we 
find that within fitting errors C ~ (aC«)i/^ Given the 
quantities modelled by the associated diffusive terms (see 
equations H31|) - (|33p ). this result is not entirely surprising. 

On the other hand. Fig. [3] reveals an important caveat 
of the closure model when applied to Rayleigh-Benard con- 
vection. The universal solution for the two horizontal stress 
components rxx{ri) and ryy{r]) can easily be deduced from 
rxx = Vyy = 0.5(r — Tzz). These horizontal stresses should 
therefore be identical to one another and have the same 
power- law dependence on as r and r^^, close to the wall 
and far from the wall. However, Fig. [3] clearly shows that 
the numerical data is at odds with the model. We attribute 
the discrepancy to the presence of large-scale coherent con- 
vective plumes in the system, which span the entire do- 
main and create strong horizontally correlated fluctuations 
as they crash against each boundaries. As a result, the fluid 




Figure 3. Comparison of the universal "convection from a wall" 
solution with numerical data for the dimensionless Reynolds 
stress components r^z, r^x and Tyy. The large symbols represent 
Tzzirj) for Ra = 2.1 X 10^ (triangles) and 2.1 X lO'' (diamonds). 
The two sets of smaller diamonds show rj:x{ri) and Tyyirf) for the 
case where Ra= 2.1 X lO'^. Note that theoretically these should 
be lying on the same curve - the difference can be attributed 
to limited statistics. In all cases Pr = 1. The dotted line shows 
the asymptotic solution r^z = rozzV^^^ using the value of Ci 
estimated by GO05, while the solid line shows a numerical inte- 
gration of the full universal profile, for our estimated parameter 
values as listed in ||46)| . (07}, and lf48] l. 



10^ 




Figure 4. Comparison of the predicted dimensionless convective 
heat flux profile /(r;) with our numerical data. The symbols have 
the same meaning as in Fig. |3] Note that the scatter for rj < 0.1 
comes from imperfect statistics very close to the wall. This plot 
was used to fit C^^ to capture the near-wall solution correctly. 
The solid line shows a numerical integration of the full universal 
profile, for our estimated parameter values as listed in II46I I. I I47II . 
and l l48t . 

in the viscous sublayer is much more strongly anisotropic 
than predicted. 

3.4-4 Far-field solution and estimation of and Ct . 

Fig. |3] compares the predicted profile for rzz iji) with data 
from our numerical simulations. The dotted line shows the 
model prediction for the solution far from the wall rzz = 
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10^ 




10"' I 

0.01 0.10 1.00 10.00 100.00 
V 

Figure 5. Comparison of the predicted dimensionless tempera- 
ture variance q with experimental and numerical data. The open 
symbols represent the results of our numerical simulations (Pr 
= 1) for Ra = 2.1 X 10^ (triangles) and 2.1 X lO'^ (diamonds). 
The plus symbols are experimental data from DuPuits, Resagk & 
Thess (2007) for Ra = 8.14x10* for air (Pr = 0.7) in a cylindrical 
box at aspect ratio 11.3. The discrepancy between the numerical 
solutions and the experimental data is attributed to the differ- 
ence between periodic side- walls and impermeable side- walls. The 
near-wall solution was used to fit Ck while the far-from-the-wall 
data was used to provide a constraint between Cq and C'7. The 
solid line show a numerical integration of the full universal profile 
as in Figs. |3] and [5] for Pr = 1. 

rzzov'^^^- Note that r^zo depends only on two numbers, the 
Prandtl number (which is known) and the model parameter 
Ci. It is reassuring to see that the value of Ci estimated by 
GO05 from wall-bounded shear flow data adequately fits the 
far-from wall solution for r^^ in this convection problem. 

The universal profiles away from the wall listed in equa- 
tion (|40|) can also be used in conjunction with numerical and 
laboratory experiments to constrain Ce and C7. These con- 
stants are unfortunately difficult to extract directly from our 
numerical simulations. The highest Rayleigh number avail- 
able (Ra — 2.1 xlO^) only has a short asymptotic (ri 3> 1) 
range, so that estimates of Ce and C7 from these datasets 
are unreliabl^U- The rms temperature data measured in var- 
ious laboratory experiments at higher Rayleigh number pro- 
vides a more adequate point of comparison. We use the rms 
temperature data of the highest aspect ratio experiments 
of DuPuits, Resagk & Thess (2007), for Ra = 8.14 x 10**. 
This dataset exhibits a significant asymptotic range, with a 
power law close to the one predicted by the closure model 
(g ~ qoV~^^^)- Fitting the data yields qo = 0.95 ± 0.05, 
which provides a first constraint between Ce and C7 (see 
Fig. [6]). Note that other datasets (from Maystrenko, Resagk 
& Thess, 2007, for example) are generally consistent with 
this estimate for qo. 

A second constraint between Ce and C7 is obtained by 
comparing the model predictions with experimental mea- 
surements of Nu(Ra). The closure model implies that Nu = 
1 -I- KKa^^^ where the constant K is a function of the 

^ This statement can be verified using a simple test problem in 
which artificial data are created using the closure model, and then 
used blindly to reconstruct Ce and C7. 




0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 

Figure 6. Calibration of the constants Cq and C7. The straight 
lines show the relationship between Cq and C7 when the con- 
stant go is equal to 0.95 (solid line), 0.9 or 1.0 (dashed lines, top 
and bottom respectively) . The curves show the value of K in the 
relationship Nu ~ 1 -|- KRa}^^, as predicted by numerical inte- 
grations of the closure model equations l|34| l for no-slip boundary 
conditions. The area marked by the intersection of the 4 dashed 
lines, and centred on the point where the two solid lines cross, 
provides estimates for C'e, and C7. 

model parameters (and the Prandtl number). The data from 
Fiinfschilling et al. (2005), Niemela & Sreenivisan (2006) and 
DuPuits, Resagk & Thess (2007) are reasonably well approx- 
imated by taking K = 0.06 ± 0.003. Variations of K with 
Prandtl number, for the range of experiments discussed, are 
within the errorbars. Given that Ci, C2, C^, Ck and C^k, are 
now known, for fixed Prandtl number, fitting K provides a 
unique relationship between Ce and C7, as seen in Fig. [H] 

By combining these two constraints, we conclude that 
a good fit to the data can be obtained with 

Ce = 1.4±0.1 , C7 = 1.4±0.1. (48) 

The values for {d} quoted in equations (|46)) . (|47)) . and ((48l) 
form from here on our selected set of parameters. These 
values are to be taken as indicative estimates, rather than 
precise calibrations. We note that the parameters derived 
do satisfy realizability (see Appendix A). The solid lines 
shown in Figs. [S] |4] and [5] are the universal boundary layer 
profiles calculated using these parameters, and are seen to 
fit all datasets (except for r^x and Vyy, as discussed above) 
satisfactorily. 

Fig. [7] compares our closure model prediction for the 
Nu(Ra) relationship, using the estimated parameters, with 
various available datasets for large aspect ratio experiments 
(F ^ 4). It also shows (as dashed lines), for comparison, 
strict upper bounds obtained by Plasting & Kerswell (2003) 
and by lerley, Kerswell & Plasting (2006) for transport by 
convection at finite and infinite Prandtl numbers respec- 
tively. It is reassuring to see that the Pr — >■ 00 prediction 
from our own closure model remains below the strict upper 
bound for the same limit. 
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Figure 7. Comparison of tlie model predictions with data for the 
Nusselt number as a function of the Rayleigh number. The square 
symbols are experimental data from Niemela &; Sreenivisan (2006) 
with Pr ~ 1 (Helium), and aspect ratio F = 4. The diamond sym- 
bols are the data from Fiinfschilling et al. (2005) with F = 6, Pr 
= 4.38 (water). The triangles are data from DuPuits et al. (2007), 
for 4 ^ F ^ 11.3, for Pr = 0.7 (air). The plus symbols are numer- 
ical data from Hartlep et al. (2007), with F = 10 and for Pr = 
0.7. Finally, the star symbols are our own numerical simulations. 
The various thin lines shows the closure model predictions for 
fiducial values of the parameters d, for Pr = 1 (solid line), Pr 
= 4.38 (dashed line) and Pr — >■ oo (dotted line). In addition, the 
two thick solid lines correspond to strict upper bound limits: the 
Nu= 1 + 0.133 Ra^/"^ line is a strict upper bound obtained by ler- 
ley, Kerswell & Plasting (2006) for Rayleigh-Benard convection 
at infinite Prandtl number, while the Nu = 1-1-0.0264 Ka}^^ line is 
a strict upper bound obtained by Plasting & Kerswell (2003) for 
Rayleigh-Benard convection at arbitrary (finite) Prandtl number. 

In conclusion, our model successfully reproduces most 
measurable features pertaining to laboratory and numerical 
experiments of Rayleigh-Benard convection, for reasonable 
values of the model parameters {Ci}. Furthermore, compar- 
ison of the estimated parameter values across a range of 
experiments in other systems (pipe flows, Couette-Taylor 
flows) shows that they are indeed of a universal nature, a 
results which can only increase confidence in our approach. 

4 HOMOGENEOUS RAYLEIGH-BENARD 
CONVECTION 

4.1 Introduction 

Another system that is of interest, and possibly more rele- 
vant to astrophysical applications, consists of an unbounded 
layer in which there is no mean flow, while the mean tem- 
perature gradient VG) is uniform and parallel to the gravita- 
tional acceleration (taken to be in the z-direction). The evo- 
lution of perturbations to this mean state can be described 
by the following set of Boussinesq equations: 

+ u ■ Vu' = —aO'gz Bz — Vtp' + vV^u , 

de' , , I de 2^/ 
at az 

V • w' = , (49) 



where all perturbations are triply periodic, as for example 

u'{x,y,z,t) = u'{x + Lx,y,z,t) 
= u'(x,y + Ly,z,t) 

= u'{x,y,z + Lz,t). (50) 

This model setup is now commonly referred to as Ho- 
mogeneous Rayleigh-Benard (HRB) convection (Borue & 
Orszag 1997; Lohse & Toschi 2003; Calzavarini et al. 2005; 
Calzavarini et al. 2006). While this system cannot be stud- 
ied using laboratory experiments, it lends itself relatively 
easily to numerical experimentation using spectral methods 
in particular. The relevant dimensionless parameters are the 
Prandtl number Pr = v/n, the Rayleigh number, now de- 
fined as 

Ra = ^ dz 5^ 

VK 

and the aspect ratio(s) F = L^^y/L^. 

The microscopic diffusivities are included in the original 
equations (|49p to regularize the system by allowing for dis- 
sipation and irreversibility. However, note that the periodic 
boundary conditions forbid the formation of boundary lay- 
ers, so it may be conjectured that the macroscopic statistical 
properties of the turbulent convection should be well defined 
and independent of v and k in the limits Ra — >■ oo (Spiegel 
1971). Furthermore, we may expect the turbulence to be 
statistically steady and homogeneous, although anisotropic. 
These properties have been argued to be more relevant to 
convection in astrophysical systems than standard Rayleigh- 
Benard convection. The HRB model may therefore provide 
a suitable local model of convection deep inside a star or 
planet. 

On dimensional grounds, the rms turbulent velocity, for 
example, must be expressible in the form 

(«5.^)'^'i./(Ra,Pr,r), (52) 

where / is a dimensionless function. According to the discus- 
sion above, / should tend to a non-zero function of F alone 
in the limit Ra oo. It is tempting to conjecture that / also 
becomes independent of F in the limit of large aspect ratio, 
F — >■ oo. This would imply that the vertical length-scale Lz 
plays a fundamental role in determining the saturation level 
of the turbulent convection, presumably by limiting the size 
of coherent structures ('eddies'). For convection deep inside 
a star or planet, it is the pressure scale-height that imposes a 
characteristic vertical scale on the turbulence (see Section[5]); 
in the local model, the vertical extent of the box plays an 
equivalent role. In practice, owing to some peculiarities of 
the HRB system discussed below, the role of the aspect ratio 
in the behaviour of the solutions is not so straightforward. 



4.2 Closure model for HRB 

4-2.1 Governing equations 

Applying our closure model to HRB, and noting that all 
statistical averages are now independent of position, we ob- 
tain the system of ODEs for the temporal evolution of the 
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second-order correlations Rij, Fi and Q: 



that there is only one positive fixed point with 



dtRxy 



Cl + C2 pl/2 p 
^ ^xy, 



dtRxz + aFxQz 



dtRyy = — -R^^^Ryy + 3^^^''^' 



dtRyz + aFyQ^ = 
dtR,, + 2aF,g, = 



Cl + C2 pl/2 p 

W -Kyz, 



3L 



dtFx + Rxz—r~ — — 7~R ^ Fx 
dz L 



P de_ Ce^i^2p 



dtFy + iiy 



Az L 

de 

Az 



d,Q + 2Fz^ = -^R'/'Q. 



Az 



L 



(53) 



where we have ignored for simplicity contributions from 
terms including Cv, Ck and Cvk which do not contribute 
to the high-Rayleigh number dynamics of HRB convection. 
Note that the resulting equation for R is 



Cl ,=,3/2 



dtR + 2aFzgz = — j^R 



(54) 



so that these equations consist of a main system 
for {R,Rzz,Fz,Q), decoupled systems for {Rxz,Fx) and 
{Ryz,Fy), and prognostic equations for Rxx,Ryy and Rxy 



4-2.2 Choice of L and consequences for the coefficients 
{C.} 

While selecting L as the distance to the wall is a natural 
choice for wall-bounded convection or shear flows, a differ- 
ent approach must be used for triply periodic flows. The 
largest eddy size in this case is limited by the horizontal 
and vertical scales in the box, so that L can be assumed to 
be proportional to min(L2:, Ly,Lz). 

It is important to note that the selection of a differ- 
ent L implies a potential rescaling of the {Ci} coefficients. 
For example, had we selected L = z/2 in the wall-bounded 
case instead of L = z, then the estimated Ci, C2, Ce and C7 
would all be half the values quoted in Section [3. 4l since these 
parameters enter the model in the combinations Ci/L, etc. 
Nevertheless, the ratios of any pairs of constants within the 
group {Cl, C2, Ce, C7} should (presumably) be preserved. 
Following these considerations, we elect to keep the esti- 
mated values of the {Ci} given in equations and (|15)) . 
and calibrate instead the value of the proportionality con- 
stant 5 in the expression L = S min(La;, Ly,Lz). 



4-2.3 High Rayleigh number HRB convection 

A search for non-trivial fixed points of the dynamical system 
(|53[) (with _R > 0) reveals they are the (positive) solutions of 
a quartic equation. In the limit of large Ra it can be shown 



Rxx — Ryy — 



C2 



Rz 



C1 + C2 
3Ci +C2'\ R 
3 ' 



R 
3' 



Cl + C2 

Rxy Rxz Ryz 

= _ Ciff/^ de 

" ^ "2L(-iV2) dl' 



F — F 

J- X — J- V 



0, 
CiR 

C7(~N^) 



de 

Az 



with 



R 



CiCe 

Note that, in this case 



Cl 3C1+C2 
Ct 3(Ci + C2) 



Q 



CeCy 



Ch 3Ch + C2_ 
Cr 3(Ci + C2) 



^2fde 

Az 



oc IVTl- 



(55) 



(56) 



(57) 



This solution represents a state of fully developed turbulent 
convection, which is statistically steady and homogeneous. 
The solution exists in the statistically ajdsymmetric sub- 



space in which Rx 



Rx 



Ry 



Fx 



Fy = 0, and is stable with respect to perturbations trans- 
verse to this subspace. It has the desired properties that 
the vertical motion is dominant {Rzz > Rxx ~ Ryy), while 
the heat flux is purely vertical and directed down the tem- 
perature gradient. Moreover, numerical integrations suggest 
that, where it exists, this state is stable and universally at- 
tracting. 

Defining the Nusselt number Nu as the ratio of the total 
to the conducted heat flux, 



Nu = 



Fz 



(58) 



we have, in the limit Ra 2> Pr 
1 



Nu 



V2C1 



Cl ^ 3Ci + C2 



X (PrRa) 



CiCs \Ct 3(Ci +C2) 
Lz 



3/2 



1/2 



(59) 



This scaling recovers the "ultimate turbulence" regime, 
where the turbulent transport properties are independent 
of microscopic diffusivities (Spiegel 1971). Defining the tur- 
bulent Reynolds number Re as Re — LF{^^^ /v, we have 



Re = 



Cl 3Ci-HC2 
CiCe \C7 3(Ci -F-C2) 



1/2 



Ra 
Pr 



1/2 



(60) 

again reproducing the standard scaling for the ultimate 
regime of convection. 



4.3 Comparison with numerical experiments 

Numerical simulations of HRB convection were first per- 
formed by Borue & Orzag (1997). More recently, Toschi & 
Lohse (2003) and Calzavarini et al. (2005) performed a range 
of Lattice-Boltzmann simulations in a cubic geometry, for 
various values of the Rayleigh and Prandtl numbers, and 
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report on the first evidence for scalings consistent witli tfie 
"ultimate regime" of convection, namely Nu oc (RaPr)^'''^ 
and Re cx (Ra/Pr)^/^. 

However, it is now recognized that the dynamics of HRB 
convection are more subtle than previously thought. As dis- 
cussed by Calzavarini et al. (2006), simulations at unit as- 
pect ratio show huge fluctuations in the instantaneous Nus- 
selt and Reynolds numbers arising from the intermittent or 
quasi-periodic (depending on Ra) exponential growth of so- 
called "elevator modes". These modes are thus named be- 
cause they are independent of z, and have the peculiar prop- 
erty of being exact nonlinear and exponentially growing so- 
lutions of the governing equations (|49p . The most unstable 
mode has a horizontal wavelength equal to the larger hori- 
zontal dimension of the box. Hence, the aspect ratio of the 
system directly influences the macroscopic solution. 

This phenomenon has a close parallel in shearing- 
box studies of the magnetorotational instability. In that 
case, forcing by a constant velocity gradient plays the role 
of the constant temperature gradient, while perturbations 
to the background fields are also assumed to be triply 
periodic. This system is unstable to equivalent "channel 
modes", exact nonlinear and exponentially growing solu- 
tions of the equations and associated periodic boundary con- 
ditions (Goodman & Xu, 1994). In this case, it is known 
that the channel modes are themselves subject to secondary 
shearing instabilities which limit their growths. However, the 
existence and growth rates of shearing instabilities depend 
sensitively on aspect ratio: they are strongly inhibited in 
systems where the streamwise direction is smaller than the 
cross-stream directions. As a result, systems with roughly 
cubic geometry are dominated by the channel modes and 
are found to have very strongly fluctuating large-scale trans- 
port properties, but for larger aspect ratio the fluctuations 
are much smaller and the channel modes are inhibited (Bodo 
et al. 2008). 

For these reasons, we performed a series of HRB simula- 
tions of various aspect ratios, in order to determine whether 
the same phenomenon occurs, and to provide a better point 
of comparison for the closure model. Appendix C provides 
a brief description of the numerical algorithm used, and the 
results are summarized in Fig. [S] We studied 5 cases, with 
= Ly and L^/L^ =1/2, 2/3, 9/10, 1/1 and 4/3. In the 
last case, the elevator modes continue growing unaffected 
by perturbations until the code fails, which seems to cor- 
roborate the premise that the secondary instabilities are in- 
hibited in wider-than-tall boxes. For F < 1, the measured 
Nusselt number eventually converges to a meaningful aver- 
age and is found to scale as predicted by the closure model, 
namely proportional to (Pr Ra)^^^F^. A good fit with the 
model predictions is found by selecting L = SLx = Lx/^/n- 
For the purpose of illustration, a snapshot of the temper- 
ature field for our largest Rayleigh number, Ra = 5 x 10® 
(with Pr = 1) and aspect ratio 1/2 is shown in Figure [S] 



4.4 The effect of rotation on homogeneous 
turbulent convection 

We now consider the effect of rotation on HRB convection, 
where the rotation axis lies at an angle 7 from the vertical 
direction: ft — (0, S7 sin7, cos 7). In this section it is more 



000 










4 1/1 ; 




1/1 




1 00 


M/iu* 1/2 






: 2/3/* 










1 







10 100 1000 10000 

(Pr Roy/' (LjLf 

Figure 8. Variation of the Nusselt number with rescaled Rayleigh 
number for Pr = 1 for homogeneous convection. The diamond 
symbols show the data from our 3D HRB numerical simulations 
for Ra = 5 X 10® and the stars for Ra = 2.16 X 10^. In all cases 
Pr = 1. The error bars show the measurement uncertainty due 
to the finite integration time of the simulation. The aspect ratio 
r = Lx/Lz of each simulation is indicated near the corresponding 
symbol. The solid line shows the asymptotic analytical solution 
II59I I. using the values of the parameters {Ci} as listed in equations 
II46I I. I I47I I. and 1 148 I I. A good fit to the data is found by choosing 
L = LxjyfK. 




Figure 9. Volume-rendered visualization of the temperature field 
for Ra = 5 X 10® and Pr = 1 for homogeneous convection in a box 
of aspect ratio 1/2. Note how, even at this high Rayleigh number, 
the size of the dominant structures is equal to the box size. 
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convenient to work with dimensionless variables so we select 
the following scalings: 



az 



gk = ggk, (61) 

where for convenience iV is defined as iV^ — — 7V^, and is 
positive when the fluid is convectively unstable. The convec- 
tive Rossby number is then defined as 



Ro = N/Q.. 



(62) 



Stationary solutions of the closure model far from onset 
of convection satisfy the following equations: 



2Ro ^{e^kiRij + ejkiRu)0.k + gtFj + gjF^ 

R 



CiR^'^R^, - C2R 



1/2 



Ri 



(63) 









7 = : 
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Figure 10. Variation of R witli Ro~^, for various values of 7, 
for values of tlie {Ci} parameters given in II46I I and II48I I. The 
Ro^^ — and Ro~^ — 00 asymptotes satisfy equations I I67II and 
II70I I respectively. 



+ 2 Ro Eijfcfijjfc 



-CaR^'^F, 



2F^ =C7R 



1/2, 



(64) 
(65) 



In the infinite Rossby number limit (equivalently in the 
non- rotating limit), the solution of these equations reduces 
to the non-dimensional form of (jSSp and (|56p . Should all of 
the quantities be expanded in terms of the inverse Rossby 
number as (for example) 



(2) 



then we find that 

7?: 2 



CiCe 



Ci 3C1 + C2 
C7 3(Ci +C2) 



+ 0(Ro"^) 



(66) 



(67) 



(and similarly for all diagonal components of R). Our ex- 
pressions for the non-diagonal terms, to first order, recover 
the equivalent of the well-known A-effect (see Riidiger, 1989) 
in the coefficient R^z'- 



Rxz 



R(0)iR 



r(0)n 



C6i?(°)(Cl +C2) 



■ sin 7 Ro 



+0(Ro" 



(68) 



The A-effect, as seen in the above equation, describes how 
rotationally constrained turbulent motions can drive dif- 
ferential rotation, through the non-diagonal component of 
the stress-tensor Rxz- As expected from dimensional analy- 
sis and geometrical arguments, its amplitude scales linearly 
with sin 7 SI. The other two components Rxy and Ryz only 
become important for more rapidly rotating systems as they 
are both 0(Ro-^). Finally, a non-negligible horizontal heat 
fiux is generated in the direction of x g, namely 



- R^xJ) + Ci + C2 Vi?(0)Ff' 

Fx = 2 i sin7Ro~^ 

l-CeR(o){Ci + C2) 

+0{Ro-^) , (69) 
although note that when applied to stellar convection zones. 



this effect is relevant only for non-axisymmetric heat trans- 
port. The "latitudinal" heat flux Fy on the other hand is of 
higher order in Ro~^. 

In the opposite limit of very low Rossby number (the 
rapidly rotating limit) an expansion in powers of Ro reveals 
that 



R 



2 cos 7 
CiCe 



Ci 3C1+C2 
Ct 3(Ci+C2) 



+ 0(Ro) 



(70) 



so that the rms velocity is reduced by a factor cos 7 com- 
pared with the non-rotating case. Note, however, how the 
expected reduction (and potential suppression) of the con- 
vective heat fiux in rapidly rotating systems where gravity is 
aligned with the rotation axis (Chandrasekhar, 1961) so that 
7 = is not captured by this closure model. This problem, 
which was identified by Miller & Garaud (2007) , can presum- 
ably be attributed to the incomplete modeling of the effects 
of the pressure-strain correlations which are known to play 
an important role in the limit of rapid rotation. It is there- 
fore likely that these effects also cause our model to yield 
inaccurate predictions for 7 7^ in the same limit. A full res- 
olution of the issue must eventually involve the derivation 
of a better closure for the pressure-strain correlation terms. 
For completeness note that in this limit the model predicts 
that a significant heat fiux is carried horizontally along By, 
with amplitude Fy — tan 7 Fz , and that 



Ry 



Ci 



C1 + C2 



sin 7 cos "/ R + 0(Ro) 



(71) 



while Rxy and R: 



are both 0(Ro). 
Fig. [To] shows the variation of the normalized as a 
function of both 7 and Ro~^, while Fig. II II shows the varia- 
tion of the normalized —Rxz/R as a function of both 7 and 
Ro~^, illustrating the dependence of the A-effect on both 
parameters as predicted by our model. 
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Figure 11. Variation of —Rxz/R with Ro~^, for various values 
of 7, for values of the {Ci} parameters given in 14611 and 14811 . 
Note that Rxz/R <^ ^ for low rotation rates, and to f2^^ for 
large rotation rates. 



4.5 Comparison with previous second-order 
models 

We now compare our findings with tlie commonly used 
model for convective stresses originally proposed by Riidiger 
& Kitchatinov (1993) and later extended by Riidiger et al. 
(2005, Ral05 hereafter). Note that the related theory of 
Kitchatinov & Riidiger (2005) relies on the presence of a 
background density stratification to explain the A-effect. As 
such it is not an appropriate point of comparison for our 
Boussinesq calculation. 

Riidiger & Kitchatinov (1993) and Ral05 assume the 
presence of a "background" turbulence caused by a given 
(unspecified) forcing mechanism, which, in the absence of 
rotation, is described by an eddy turnover time r, a mix- 
ing length / and a turbulent diffusivity Vt = I'^/t. This 
background turbulence also may also have some degree of 
anisotropy, controlled by the parameter a defined in our no- 
tation as 



r{0) , r(0) 



2R 



(0) 



R 



(0) 



(72) 



where the superscript (0) denotes turbulent quantities of 
the non-rotating system. Note how a — for isotropic tur- 
bulence. 

Ral05 show how the presence of rotation (where the 
rotation axis lies at an angle 7 from the vertical) modifies 
the background turbulence, an effect which gives rise to non- 
diagonal components in the stress tensor. They argue that 
this phenomenon is controlled by the Coriolis number il* 
defined as 



(73) 



Their eddy turnover time r is naturally related to L/y/Rf^ 
in our notation, so that, for the purpose of comparison we 
have 



Q,* oc 



where the proportionality constant is of order unity. 



(74) 



In the limit of slow rotation, Ral05 predict a A— effect 
through the following term: 

2a . QL 



■ sm 7- 



-.R 



(0) 



(75) 



where the proportionality constant is the same as in equa- 
tion (|74l) . Meanwhile, our model when written in dimen- 
sional form and recast in terms of the anisotropy factor a 
yields 



Rx 



Ci 



3(Ci -f C2) 



-C6i?(0)(Ci +C2) 3Ci +C2 



sm7 ■ 



R 



(0) 



(76) 

where R^'^ is a dimensionless constant which depends only 
on the model parameters, and is given by equation (|67|) with 
Ro~^ = 0. In the same slow-rotation limit, the other off- 
diagonal components of the stress tensor are 0(Ro~^) or 
higher order in both our and their models. 

Overall, the two formalisms agree on the dependence 
of the stresses on the rotation rate and on latitude, as ex- 
pected on dimensional and geometrical grounds. In addi- 
tion, both models explicitly demonstrate the importance of 
the anisotropy of the background non-rotating turbulence 
in controlling the amplitude of the A-effect. However, the 
dependence of Rxz on the anisotropy factor a superficially 
appears to be different in the two theories. We interpret 
this in two ways. First, note that the anisotropy factor a 
is a "free" parameter in the works of Ral05. In our model 
by contrast, there is no freedom in independently specify- 
ing a since it is a solution of the model once the system is 
specified (e.g. shearing flow, convective flow) and depends 
on the {Ci} parameters. In the HRB system for example 
a = -6Ci/(3Ci + C2). 

Secondly, Rxz is directly proportional to a in the model 
of Ral05 while our model reveals an additional contribution 
to the A-effect arising from the background turbulent heat 
flux (see equation (|68|1 for a more explicit expression). This 
contribution is missing from the model of Ral05 which does 
not take into account the heat equation. As a result, one may 
superficially conclude that the A-effect could exist even for 
isotropic background turbulent convection. In practice, it is 
difficult to conceive of a naturally occurring isotropic turbu- 
lent system which has a non-zero vertical heat flux, so the 
term is in fact also indirectly related to the anisotropy 
of the system, although perhaps not exactly in the same way. 

Finally, we emphasize that in the limit of rapid rota- 
tion, neither theory is expected to be accurate because of 
the extreme induced anisotropy of the rotating turbulent 
motions. Nevertheless it is interesting to note that the pre- 
dicted dependence of the stresses on the rotation rates now 
no longer agree with one another. We find that Ryz tends 
to a constant independent of rotation rate while Ral05 find 
that Ryz oc Ro. For the other off-diagonal components Rxz 
and Rxy we find a dependence on Ro, while they predict a 
dependence on Ro'^. 

We conclude this section by emphasizing the success of 
our closure model in reproducing numerical experiments of 
HRB convection at various aspect ratios and Rayleigh num- 
bers. Furthermore our model predictions are exactly pro- 
portional to those of Ral05 (with a proportionality constant 
of order unity) for convection in a slowly rotating system. 
Hence we expect to recover many of the results and suc- 
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cesses of these authors in modehng differential rotation in 
stars, albeit with an extended model which self-consistently 
includes heat transport in addition to angular momentum 
transport. In preparation of this future modeling endeav- 
our, we finally turn to the next natural step of this work, 
namely the extension of the model to the anelastic and fully 
compressible equations. 



5 THE ANELASTIC SYSTEM AND 
COMPRESSIBLE FLOWS 

So far we have worked within the Boussinesq approxima- 
tion, which is applicable only to a shallow layer of fluid 
whose depth is much less than the density scaleheight. In 
order to apply our model to stars we must first adapt it to 
the anelastic approximation (Ogura & Phillips 1962; Gough 
1969), which is relevant to subsonic convection in a deep 
layer. 

Here we follow the more standard derivation of the 
anelastic approximation where the reference state is taken 
to be an adiabatically stratified fiuid in hydrostatic equi- 
librium. The reference density po{r) and temperature To{r) 
may vary substantially, while the specific entropy sq is uni- 
form. In place of equations (P))- (|ll|l we have 



[dt + Ujdj)ui = — (s — 3o)diTo — diTp + ■ 



{dt + Uidi){s — so) = 



(77) 
(78) 



(79) 



where the dots represent terms due to viscosity (in the equa- 
tion of motion) and thermal conduction (in the thermal en- 
ergy equation), while V" is, again, a modified pressure. Vis- 
cous dissipation can also be included in the thermal energy 
equation, although it is usually omitted in the Boussinesq 
approximation. A derivation of these equations, omitting dif- 
fusive effects, is given in Appendix B. 

The anelastic system is formally very similar to the 
Boussinesq system except for the variable density of the ref- 
erence state. However, the entropy perturbation and back- 
ground temperature gradient respectively play the roles 
taken by the temperature perturbation and agi in the 
Boussinesq approximation. A very similar analysis to that 
carried out for the Boussinesq system leads to equations for 
Rij, Fi and Q of the form 

{dt + Ukdk)Rij + RikdkUj -\- RjkdkUi + RijdkUk 

+FidjTo + F,d,To = , (80) 



{dt + Ujdj)Fi + RijdjS + FjdjUi + FidjUj + QdiTo 



{dt + Uidi)Q -f 2FidiS 4- QdiUi 



(81) 
(82) 



where the dots represent terms that require a closure model. 
In the anelastic system the relevant definitions of the 
Reynolds stress Rij, flux Fi and variance Q are 



Rij = {pv,u[u'j), 



Fj = {pou.s ), Q = {pos ). (83) 

Note that Rij now has the correct dimensions for a stress 
tensor, and that Fi is really an entropy flux density. Some 



additional linear terms arise in the anelastic system because 
diUi / 0. 

We apply the same closure model as for the Boussi- 
nesq system, except that the relaxation timescale which was 
proportional to L/R^^'^ is now proportional to L/{R/ poY^'^ 
because of the redefinition of Fiij : 

{dt -I- Ukdk)Rij + RikdkUj -f RjkdkUi + RijdkUk 



+Fd,To + F,(9.To = ( — 



L \po 



1/2 



Ri 



Ca fR 
Po 



1/2 



(84) 



{dt + Ujdj)Fi -f RijdjS + FjdjUi + FidjUj + QdiTo 

1/2 _ 

F,, (85) 



Ce fR_ 
L \po 



{dt + uA)Q + 2FAs + Qd^U, = -^( — 

L \po 



1/2 



(86) 



We do not include any of the terms proportional to or k 
here because we consider the high-Rayleigh number limit in 
the absence of rigid boundaries only. 

The question arises as to how the length-scale L should 
be identified for anelastic convection in a deep layer. It 
should probably related to the pressure scaleheight or den- 
sity scaleheight, as in the stellar mixing-length theory. In- 
deed, numerical simulations of convection in spherical shells 
with a substantial density variation indicate that the convec- 
tive cells are much smaller near the outer surface where the 
scaleheight is small; nevertheless, there may be situations in 
which convective plumes can span several scaleheights. 

Equations (|84p - (|86p can then be combined with equa- 
tions for the mean variables in the form 



di{poUi) = 0, 



po{dt + Ujdj)ui = ~{s ~ so)a,To - podii) - 



(87) 



djRij, 



poTo{dt + uA){s~so) = ^(—] ^~Tod^F. (89) 

L \po J 2 

In the last equation we have included the turbulent viscous 
heating. 

These equations could be applied to studying convec- 
tion and meridional circulation in rotating stars. The solu- 
tion can be assumed to be axisymmetric and independent 
of time, although for practical purposes it may be easier to 
evolve the equations forwards in time until a steady state is 
reached (if it is) rather than directly seeking such a solution. 

In the absence of rotation the problem becomes spher- 
ically symmetric, the mean flow disappears, the stress be- 
comes diagonal (although anisotropic) and we obtain the 
local algebraic system 

2F.drTo = -^1+^ ( I) R^^ + %.( ^] R, (90) 
L \po J 6L \po J 



1/2 



R 



2f;arTo 



CTi 
' L 



R^ 



R, 



(91) 
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RrrdrS + QdrTo = — ^ ( — 



1/2 



Cr f R 



1/2 



(92) 



(93) 



The solution is, by direct analogy with equations (|55|l - (|56p . 



Reo ~ R, 



V + Ca y 3 ' 

Ci + CaJ 3' 



_ C^iR/pof^ 

— ;7T7 — JTol — PoOrS, 



2L[-m) 



Q 



CiR 



{drS) , 



with 



R = 



CiCe 



Ch^ 3C1 + C2 



C7 3(Ci+C2) 



(94) 



poi'(-iV^), (95) 



where now 



= [drTo)drS. 



(96) 



In this situation the entropy gradient s is not known in ad- 
vance. However, to balance the thermal energy equation, 
diFi = 0, which implies that r^F,. is a constant, determined 
by the luminosity generated by the core of the star. (This 
conclusion is modified if the radiative flux or any sources 
of energy such as turbulent viscous dissipation make an im- 
portant contribution to the thermal energy equation.) Then 
the above equations can be solved algebraically to find drS, 
R, etc., at each radius, assuming that a prescription for L is 
given. The result is equivalent to a version of mixing-length 
theory. 

Rotation couples radial and latitudinal transport of 
heat and momentum and induces large-scale entropy gra- 
dients and mean fiows. However, if we assume that their 
effects can be ignored in the overall turbulent dynamics con- 
trolling the properties of the stresses, then the local A— effect 
is easily recovered as an anelastic version of equation (|68p . 
As before, the only differences with the Boussinesq case is 
that (i) the two terms containing T?'"' , which have their ori- 
gin in the eddy turnover time, should be replaced by R}-^^ / pa 
and (ii) in expressing H68p in dimensional form (see equation 
l|6ip ). one must also replace iV^ by {drTo)drS and AQ/dz by 
podrS, as seen above. The resulting expression then directly 
links the turbulent transport of angular momentum and of 
heat to one another. Since heat transport in this model is 
very similar to mixing-length theory, our formalism now pro- 
vides a simple framework in which to combine models of 
stellar structure with models of internal stellar dynamics. 
Note that in practice mean flows and especially latitudinal 
entropy gradients could play a role in the global dynamics 
of the system. The whole model should therefore be solved 
self-consistently and globally instead of using (|68p . This can 
only be done numerically and is deferred to a subsequent 
paper. 

It is also possible to 'import' the model of anelastic con- 
vection into the full set of equations governing the motion 
of a compressible fluid. The idea here is that, while the con- 
vection might be assumed to be subsonic and to obey the 



anelastic approximation, the mean flow need not obey these 
constraints. An example is convection in an accretion disc, 
where the accretion flow, although slow, cannot be treated 
in the anelastic approximation with a reference density pro- 
file. Omitting now the bars on all quantities, and neglecting 
self-gravitation (although it can easily be restored), we pro- 
pose a system of equations consisting of the equation of mass 
conservation, 

^tp + ^^{pu^)=Q, (97) 
the equation of motion, 

p{dt + Ujdj)ui = —pdi^ — dip — djRij, (98) 
and the thermal energy equation, 

Ci fRV^^ R 



pT{dt + uA)s = -J- j - - Td,F„ (99) 

together with the equations of the closure model, 

(9t -I- ukdk)Rij + RikdkUj + RjkdkUi + RijdkUk 



+F,d,T + F.d^T ^ ( - 

P 



1/2 



Ri 



1/2 



(dt + Ujdj)Fi + RijdjS + FjdjUi + FidjUj + QdiT 



(100) 



(101) 



(dt + uA)Q + 2FAs + Qd.,u, = - ^ ' Q- (102) 

The total energy is then exactly conserved in the form 
dt[p{\u +^ + e) + \R] 

+ft [p{\u $ h)u^ + \Ru^ + R^jUj + TF^] = 0, 

(103) 

where e and h are the specific internal energy and the spe- 
cific enthalpy, respectively, and the gravitational potential 
$ is assumed to be independent of time. The existence of 
this conservation law implies a certain self-consistency in 
the equations of the model. The terms that were added in 
passing to the compressible model are required to have the 
form that they do in order that energy be conserved. We 
note again that Fi is really an entropy fiux density, and that 
TFi is the corresponding energy flux density. 

The physical content of this model is that the turbu- 
lent convecting fluid behaves similarly to a complex, non- 
Newtonian material in which there is a dynamical consti- 
tutive equation that relates the stress tensor to the defor- 
mation history of the fluid. The above equation for dtRij 
(along with those for dtFi and dtQ) plays this role. 



6 CONCLUSIONS AND FUTURE PROSPECTS 

We have laid out the foundations of a new second-order clo- 
sure model for the dynamics of turbulent convection, with 
future applications to stellar convective regions in mind. 
This model is a direct extension of the work of Ogilvie (2003) 
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and GO05, and has similar properties. The proposed clo- 
sure has a straightforward physical interpretation, and well- 
understood limitations. 

Comparison with laboratory and rmmcrical experiments 
reveals good overall agreement of the model predictions with 
known properties of rotating shear flows (GO05) and high 
Rayleigh-number rotating convection (this work). In partic- 
ular, our model naturally reproduces the standard scaling 
relationships between the Rayleigh and Nusselt numbers for 
Rayleigh-Benard convection and for Homogeneous Rayleigh- 
Benard convection, and contains the well-known A-effect de- 
scribing angular momentum transport in a rotating turbu- 
lent fluid. 

When extended to the anclastic (or fully compressible) 
case, our formalism can be applied to study convection in 
stellar interiors. Note that the effects of Maxwell stresses can 
also straightforwardly be included following Ogilvie (2003) 
if needed. We show that the model naturally reduces to 
a version of mixing-length theory when applied in a one- 
dimensional framework. In the presence of rotation it be- 
comes a powerful tool to study within a single framework the 
multi-dimensional balance involving large-scale mean quan- 
tities such as the entropy profile, the meridional circulation 
and the differential rotation. Future work applying our clo- 
sure model in a spherical shell geometry will help understand 
some of the trends seen in the increasingly large number of 
available observations of stellar differential rotation. 
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APPENDIX A: REALIZABILITY 

In this Appendix we show that, provided that the condition 

2C6 - Ct - Ci - C2 ^ (Al) 

is satisfied, the evolutionary model used in this paper en- 
sures that the Reynolds stress Rij remains positive semi- 
definite and the entropy variance Q remains positive at all 
points in the flow if they are so initially. This will ensure that 
the quantities Rij , Fi and Q predicted by the model can be 
realized by genuine velocity and entropy perturbations. Not 
least, it will ensure that the turbulent kinetic energy remains 
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non- negative. We work here with the anelastic version of the 
closure model, equations H84p - (|86p , although a similar argu- 
ment applies to the Boussinesq system in the high-Rayleigh 
number limit in which the terms proportional to or «; are 
omitted. 

Let Ai be any vector with appropriate dimensions, and 
consider the tensor (within the anelastic system) 



Sij = (^po{u'i + Ais'){u'j + Ajs')^ 



(A2) 



at any point in the flow. The associated quadratic form is 

S = S.jX.Xj = (po [{u + As') ■X]^y (A3) 

Evidently 5^0 for all vectors Xi, and therefore Sij must be 
a positive semi-deflnite tensor, for any choice of Ai, at every 
point in the flow. Allowing the vector Ai to vary provides 
us with a family of realizability conditions. 

On the other hand, Sij can be expanded as 

Sij = Rij + FiAj +FjAi + QAiAj , (A4) 

and therefore 

S = R^jXaj+2(F-X){A-X) + Q{A-Xf 



{R^,-Q-^FF,)X,Xj+' 



[{F + QA) ■ X] 



(A5) 



Provided that Q > 0, a necessary and sufflcient condition 
for S to be non-negative for all choices of Ai and Xi is that 
the tensor 



Tij — Rij Q Fi Fj 



(A6) 



be positive semi-deflnite. If this condition is satisfled then 
TijXiXj ^ for all Xi and therefore 5^0 for all Xi 
and Ai. On the other hand, if a vector Xi exists such that 
TijXiXj < 0, then 5 < for this choice of Xi if we set 
A, = ^F/Q. 

We therefore aim to show that, provided that the con- 
dition (|A1|) is satisfled, the model ensures that Tij remains 
positive semi-deflnite at all points in the flow if, in the initial 
state, Tij is positive semi-deflnite and Q > at all points. 

We apply a reduction ad absurdum. If Tj has a negative 
eigenvalue at any event, then the quadratic form 

r = T.X^Xj (A7) 

is negative at that event, for some choice of the vector Xi. 
Without loss of generality, let Xi be a differentiable vector 
fleld advected according to the time-reversible equation 



DX.-Xjdii 



0, 



(A8) 



and such that T < at the event in question. Here D = dt + 
Uidi is the Lagrangian derivative following the mean flow. 
Retracing the the value of T to the initial state, following 
the mean flow in reverse, we deduce that T must have passed 
through zero with DT < 0. However, using equations (|84[) " 
(l86l) we flnd 



DT = -Td,u, + 2Q'^{F ■ X)XmjdjS 

-(Ci + OL-' ( ^\ T + C^L-' ( ^ ^X' 
\PoJ \poJ 3 

-^(2(76 - Ct - Ci - C2)L"M ^ j Q-^{F-Xf. 



When T = 0, Xi is a null eigenvector of Tj and therefore 
DT ^ 0, with equality only when there is no turbulence. 
Therefore Tij cannot in fact develop a negative eigenvalue. 



APPENDIX B: NUMERICAL ALGORITHM 
FOR RAYLEIGH-BENARD CONVECTION 

The equations governing Rayleigh-Benard-convection, (9- 
11), can be written in the standard non-dimensional form 

dtVL - PrV^u = 



dtT 



-VV" + Pr RaTe^ 
-V • (uT) , 



Vu 







(V X u) X u ,(B1) 
(B2) 
(B3) 



where is the vertical unit vector and non-dimensional 
quantities are marked by a hat. Here the layer height h is 
used as the unit length, / k is the unit time and the im- 
posed temperature difference between the two plates AT 
serves as the temperature scale. 

We decompose the velocity field u = {ux,Uy,Uz) into 
toroidal, poloidal and mean parts, 

u = V X (ee2)4-V X V X {f e^) + {u^)he^ + {uy)hey , (B4) 

where {...)h denotes a horizontal average. Note that (|B4|I au- 
tomatically satisfies (|B3|I . Equations for the scalar functions 
e and / and for the horizontally averaged velocities {ux)h 
and {uy)h can be derived by applying the operators ■ Vx 
and Gz ■ V X Vx to the momentum equation (|Bip and by 
averaging this equation horizontally, leading to 



(9t - PrV')Vlfe 
{dt - PrV')VlfV'/ 

{dt ~ Pr dl){ux)h 
{dt~VTdl){uy)u 



■ V X N , 
-e^ • V X V X N 
PrRaVlfT , 
-e^ • (N)h , 
-e« • (N)h , 



(B5) 

(B6) 
(B7) 
(B8) 



where V|/ — ~\- dy is the horizontal Laplacian and where 
the vector quantity N is defined as 



N = (V X u) X u . 



(B9) 



Either stress-free or no-slip boundary conditions may be ap- 
plied for u, which translates into the boundary conditions 

dzC — f — dlf = dz{ux) = dz{uy) = (stress free) 



e = dzf = dlf = {ux) = {uy) = (no slip) 



(BIO) 



(A9) 



at z = and z = 1 for the toroidal and poloidal scalars e 
and / and for the mean velocities (ux) and (uy). 

A pseudo-spectral algorithm is used to solve the govern- 
ing equations in the formulation HB2IB5IB8|) . The fields e, / 
and T are expanded as Fourier series in the horizontal direc- 
tion. In the vertical direction, a Chebyshev expansion on a 
Gauss-Lobatto grid is used for all unknowns. Fast transform 
algorithms can then be applied to switch between physical 
and transform space. A semi-implicit time-stepping scheme 
is employed for the temporal discretization, where all linear 
terms are treated implicitly by a second order Backward- 
Differencing (BDF2) scheme, while a second order Adams- 
Bashforth (AB2) scheme is applied to the nonlinear terms. 

Most of the computation is carried out in spectral space, 
although the nonlinear terms are evaluated in physical space. 
The usual 3/2-rule is applied to avoid aliasing errors in the 
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horizontal directions, whereas no de-ahasing procedure is 
employed along the vertical coordinate. We use the Cheby- 
shev tau method (Peyret, 2002) to solve the ODEs arising 
from the implicit part of the time-stepping scheme. This 
method has the advantage of yielding linear systems which 
can be manipulated into sparse form, thus keeping the mem- 
ory requirements at a manageable level. The code is paral- 
lelized using transpose-based parallel FFTs, see Stellmach 
& Hansen (2008) for details. 



APPENDIX C: NUMERICAL ALGORITHM 
FOR HOMOGENEOUS RAYLEIGH-BENARD 
CONVECTION 

A spectral algorithm using Fourier expansions in all three 
spatial directions is used to solve the governing equations 
(49) in the homogeneous Rayleigh-Benard case. The primi- 
tive variables u' ,T' ,ip' are used, with the pressure perturba- 
tion i/)' being calculated in the same way as in the classical 
Patterson-Orzag Algorithm (Canuto et al. 2007) which is 
widely used in simulations of homogeneous, isotropic tur- 
bulence. Non-linear products are evaluated on a grid in 
physical space and aliasing errors are avoided by using the 
3/2-rule. The equations are advanced in time by a semi- 
implicit multistep method in which all diffusive terms are 
treated implicitly by a third order Backward-Differencing 
(BDF3) algorithm, while a third-order Adams-Bashforth 
(AB3) scheme is applied to the non-linear terms. This time 
stepping method offers a relatively large stability domain 
that includes a part of the imaginary axis at a compara- 
tively low cost (Peyret, 2002). As a starting scheme, we use 
a second-order Runge-Kutta method. A parallelization ap- 
proach similar to the one employed in our Rayleigh-Benard 
convection code is used, see Stellmach & Hansen (2008) for 
details. 



APPENDIX D: DERIVATION OF THE 
ANELASTIC SYSTEM 

The equations governing the motion of an ideal, compress- 
ible fluid are the equation of mass conservation, 

dtp + d^{pu,) = 0, (Df) 

the equation of motion, 

p{dt + Ujdj)ui = —pdi^ — dip, (D2) 

the thermal energy equation, 

pT{dt + uA)s = 0, (D3) 

and Poisson's equation, 

dr,^ = A-nGp. (D4) 

By introducing the speciflc enthalpy h, the equation of mo- 
tion can be rewritten in the form 

{dt + ujdj)u, = -94$ + h)+ Td,s. (D5) 

We adopt a system of units in which the pressure scale- 
height and the sound speed are of order unity. We intro- 
duce a small parameter e such that the (imaginary) Brunt- 
Vaisala frequency is 0(e) when expressed in these units. We 



then pose the asymptotic expansions 



o 


Oa(v\ -\~ Oo(v 




n = 


€Ui{r,T) + 0{e 




$ = 


<I>o(t') + e^$2(r' 




h = 


/io(r) + e^/i2(r, 


r) + 0{e'), 


T = 


To{r) + e^T2{r, 




s = 


So + f?S2(r, t) - 


^o(.*). 



where r = is a slow time variable. Note that the reference 
state is adiabatically stratified, and therefore so is indepen- 
dent of r. The equation of motion at leading order [0(1)] 
requires hydrostatic equilibrium in the reference state, 

= -a($o + /ia), (D7) 

while Poisson's equation at leading order [0(1)] is 

du^o = 47rGpo. (D8) 

At O(e^) the equation of motion gives 

{dr + Uljdj)uu = -a($2 + h2) + Tod^S2. (D9) 

This can also be written in the form 

{Ot + uijdj)uu = ~S2diTo — dit/j, (DIO) 

where t/) = <I>2 + ^2 — roS2 is a modifled pressure variable. 
The equation of mass conservation at leading order [0(e)] is 

di{pouu) = 0, (Dll) 

and the thermal energy equation at leading order [0(e'')] is 

poToid^ + ui,d,)s2^0. (D12) 

Poisson's equation at O(e^) is not required, and the de- 
partures from the reference state are not affected by self- 
gravitation at leading order. When the asymptotic scalings 
are removed, and allowance is made for diffusive effects, 
equations (|77fl - (|79|) are obtained. 
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